﻿/*
 * This file is part of MonoStrategy.
 *
 * Copyright (C) 2010-2011 Christoph Husse
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Affero General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Affero General Public License for more details.
 *
 *  You should have received a copy of the GNU Affero General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Authors: 
 *      # Christoph Husse
 * 
 * Also checkout our homepage: http://monostrategy.codeplex.com/
 */
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace MonoStrategy
{
    internal class PerlinNoise
    {

        internal void FillNoiseArray2D(
           int inOctaves,
           double inFrequency,
           double inAmplitude,
           int inSeed,
           int inWidth,
           int inHeight,
           double[] outArray)
        {
            PerlinNoise noise = new PerlinNoise(inOctaves, inFrequency, inAmplitude, inSeed);
            double[] vec = new double[2];

            for (int y = 0, offset = 0; y < inHeight; y++)
            {
                for (int x = 0; x < inWidth; x++)
                {
                    vec[0] = x;
                    vec[1] = y;

                    outArray[offset++] = noise.Perlin_noise_2D(vec);
                }
            }
        }

         internal void FillNoiseArray3D(
             int inOctaves,
             double inFrequency,
             double inAmplitude,
             int inSeed,
             int inWidth,
             int inHeight,
             int inDepth,
             double[] outArray)
         {
             PerlinNoise noise = new PerlinNoise(inOctaves, inFrequency, inAmplitude, inSeed);
             double[] vec = new double[3];

             for (int z = 0; z < inDepth; z++)
             {
                 for (int y = 0, offset = 0; y < inHeight; y++)
                 {
                     for (int x = 0; x < inWidth; x++)
                     {
                         vec[0] = x;
                         vec[1] = y;
                         vec[2] = z;

                         outArray[offset++] = noise.Perlin_noise_3D(vec);
                     }
                 }
             }
         }


        private const int SAMPLE_SIZE = 1024;
        private const int B = SAMPLE_SIZE;
        private const int BM = (SAMPLE_SIZE - 1);
        private const int N = 0x1000;
        private const int NP = 12;
        private const int NM = 0xfff;

        private int mOctaves;
        private double mFrequency;
        private double mAmplitude;
        private int mSeed;
        private CrossRandom mRngCtx;
        private int[] p = new int[SAMPLE_SIZE + SAMPLE_SIZE + 2];
        private double[][] g3 = new double[SAMPLE_SIZE + SAMPLE_SIZE + 2][];
        private double[][] g2 = new double[SAMPLE_SIZE + SAMPLE_SIZE + 2][];
        private double[] g1 = new double[SAMPLE_SIZE + SAMPLE_SIZE + 2];

        private static double s_curve(double t) { return (t * t * (3.0f - 2.0f * t)); }
        private static double lerp(double t, double a, double b) { return (a + t * (b - a)); }

        double Noise1(double arg)
        {
            int bx0, bx1;
            double rx0, rx1, sx, t, u, v;
            double[] vec = new double[1];

            vec[0] = arg;

            t = vec[0] + N;
            bx0 = ((int)t) & BM;
            bx1 = (bx0 + 1) & BM;
            rx0 = t - (int)t;
            rx1 = rx0 - 1.0f;

            sx = s_curve(rx0);

            u = rx0 * g1[p[bx0]];
            v = rx1 * g1[p[bx1]];

            return lerp(sx, u, v);
        }

        double Noise2(double[] vec)
        {
            int bx0, bx1, by0, by1, b00, b10, b01, b11;
            double rx0, rx1, ry0, ry1, sx, sy, a, b, t, u, v;
            double[] q;
            int i, j;

            t = vec[0] + N;
            bx0 = ((int)t) & BM;
            bx1 = (bx0 + 1) & BM;
            rx0 = t - (int)t;
            rx1 = rx0 - 1.0f;

            t = vec[1] + N;
            by0 = ((int)t) & BM;
            by1 = (by0 + 1) & BM;
            ry0 = t - (int)t;
            ry1 = ry0 - 1.0f;

            i = p[bx0];
            j = p[bx1];

            b00 = p[i + by0];
            b10 = p[j + by0];
            b01 = p[i + by1];
            b11 = p[j + by1];

            sx = s_curve(rx0);
            sy = s_curve(ry0);

            q = g2[b00];
            u = (rx0 * q[0] + ry0 * q[1]);
            q = g2[b10];
            v = (rx1 * q[0] + ry0 * q[1]);
            a = lerp(sx, u, v);

            q = g2[b01];
            u = (rx0 * q[0] + ry1 * q[1]);
            q = g2[b11];
            v = (rx1 * q[0] + ry1 * q[1]);
            b = lerp(sx, u, v);

            return lerp(sy, a, b);
        }

        double Noise3(double[] vec)
        {
            int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
            double rx0, rx1, ry0, ry1, rz0, rz1, sy, sz, a, b, c, d, t, u, v;
            double[] q;
            int i, j;

            t = vec[0] + N;
            bx0 = ((int)t) & BM;
            bx1 = (bx0 + 1) & BM;
            rx0 = t - (int)t;
            rx1 = rx0 - 1.0f;

            t = vec[1] + N;
            by0 = ((int)t) & BM;
            by1 = (by0 + 1) & BM;
            ry0 = t - (int)t;
            ry1 = ry0 - 1.0f;

            t = vec[2] + N;
            bz0 = ((int)t) & BM;
            bz1 = (bz0 + 1) & BM;
            rz0 = t - (int)t;
            rz1 = rz0 - 1.0f;

            i = p[bx0];
            j = p[bx1];

            b00 = p[i + by0];
            b10 = p[j + by0];
            b01 = p[i + by1];
            b11 = p[j + by1];

            t = s_curve(rx0);
            sy = s_curve(ry0);
            sz = s_curve(rz0);

            q = g3[b00 + bz0]; u = (rx0 * q[0] + ry0 * q[1] + rz0 * q[2]);
            q = g3[b10 + bz0]; v = (rx1 * q[0] + ry0 * q[1] + rz0 * q[2]);
            a = lerp(t, u, v);

            q = g3[b01 + bz0]; u = (rx0 * q[0] + ry1 * q[1] + rz0 * q[2]);
            q = g3[b11 + bz0]; v = (rx1 * q[0] + ry1 * q[1] + rz0 * q[2]);
            b = lerp(t, u, v);

            c = lerp(sy, a, b);

            q = g3[b00 + bz1]; u = (rx0 * q[0] + ry0 * q[1] + rz1 * q[2]);
            q = g3[b10 + bz1]; v = (rx1 * q[0] + ry0 * q[1] + rz1 * q[2]);
            a = lerp(t, u, v);

            q = g3[b01 + bz1]; u = (rx0 * q[0] + ry1 * q[1] + rz1 * q[2]);
            q = g3[b11 + bz1]; v = (rx1 * q[0] + ry1 * q[1] + rz1 * q[2]);
            b = lerp(t, u, v);

            d = lerp(sy, a, b);

            return lerp(sz, c, d);
        }

        void Normalize2(double[] v)
        {
            double s;

            s = (double)Math.Sqrt(v[0] * v[0] + v[1] * v[1]);
            s = 1.0f / s;
            v[0] = v[0] * s;
            v[1] = v[1] * s;
        }

        void Normalize3(double[] v)
        {
            double s;

            s = (double)Math.Sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
            s = 1.0f / s;

            v[0] = v[0] * s;
            v[1] = v[1] * s;
            v[2] = v[2] * s;
        }

        internal PerlinNoise(int octaves, double freq, double amp, int seed)
        {
            int i, j, k;

            mOctaves = octaves;
            mFrequency = freq;
            mAmplitude = amp;
            mSeed = seed;
            mRngCtx = new CrossRandom(mSeed);

            for (i = 0; i < g3.Length; i++)
            {
                g3[i] = new double[3];
                g2[i] = new double[2];
            }

            for (i = 0; i < B; i++)
            {
                p[i] = i;
                g1[i] = (double)((mRngCtx.Next() % (B + B)) - B) / B;
                for (j = 0; j < 2; j++)
                    g2[i][j] = (double)((mRngCtx.Next() % (B + B)) - B) / B;
                Normalize2(g2[i]);
                for (j = 0; j < 3; j++)
                    g3[i][j] = (double)((mRngCtx.Next() % (B + B)) - B) / B;
                Normalize3(g3[i]);
            }

            while (--i > 0)
            {
                k = p[i];
                p[i] = p[j = mRngCtx.Next() % B];
                p[j] = k;
            }

            for (i = 0; i < B + 2; i++)
            {
                p[B + i] = p[i];
                g1[B + i] = g1[i];
                for (j = 0; j < 2; j++)
                    g2[B + i][j] = g2[i][j];
                for (j = 0; j < 3; j++)
                    g3[B + i][j] = g3[i][j];
            }

        }

        internal double Perlin_noise_3D(double[] vec)
        {
            int terms = mOctaves;
            double freq = mFrequency;
            double result = 0.0f;
            double amp = mAmplitude;

            vec[0] *= mFrequency;
            vec[1] *= mFrequency;
            vec[2] *= mFrequency;

            for (int i = 0; i < terms; i++)
            {
                result += Noise3(vec) * amp;
                vec[0] *= 2.0f;
                vec[1] *= 2.0f;
                vec[2] *= 2.0f;
                amp *= 0.5f;
            }


            return result;
        }

        internal double Perlin_noise_2D(double x, double y)
        {
            double[] vec = new double[2];

            vec[0] = x; vec[1] = y;

            return Perlin_noise_2D(vec);
        }

        internal double Perlin_noise_2D(double[] vec)
        {
            int terms = mOctaves;
            double freq = mFrequency;
            double result = 0.0f;
            double amp = mAmplitude;

            vec[0] *= mFrequency;
            vec[1] *= mFrequency;

            for (int i = 0; i < terms; i++)
            {
                result += Noise2(vec) * amp;
                vec[0] *= 2.0f;
                vec[1] *= 2.0f;
                amp *= 0.5f;
            }

            return result;
        }
    }
}
